function [err]=FP_IO(DFILE, sf, plot_flag)
% Analyze Field Potential time course

global VOLTAGE CURRENT
global CONTROL
%

err = 0;

[DFILE] = synch_file(DFILE, sf);
if(isempty(DFILE) | isempty(VOLTAGE))
   fprintf('? could not get dfile');
   return;
end;

QueMessage('FP_IO - Starting');
dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 0;	% 0 turns off detailed spike analysis; 1 turns it on 

update = 0; % flag
if(~check_field(CONTROL(sf),'stim_time')) % doesn't exist - then create it.
   CONTROL(sf).stim_time=[5 55];
   QueMessage('Updated Control with stim_time');
   update = 1;
end;
if(~check_field(CONTROL(sf),'psp_time'))
   CONTROL(sf).psp_time= [9 11 31, 9 51 71];
   QueMessage('Updated Control with psp_time');
   update=1;
end
if(~check_field(CONTROL(sf),'psp_reclist'))
   CONTROL(sf).psp_reclist= '';
   QueMessage('Updated Control with psp_reclist');
   update=1;
end
if(~check_field(CONTROL(sf),'sign_list')) % sign of what we detect (-1 for min, 1 for max)
   CONTROL(sf).sign_list= '';
   QueMessage('Updated Control with sign_list');
   update=1;
end

if(update==1)
   QueMessage('Updated Control with psp params for IO');
   return;
end
% abstract general information
QueMessage('Preparing for analysis');
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate*DFILE.nr_channel/1000;

[records,pts]=size(VOLTAGE);
% always compute the time of each stimulus, in seconds from the start of the data
a = DFILE.ztime;
TM = (a - a(1))/1000; % time at which data was collected... gathered from ztime...
IO = seq_parse(CONTROL(sf).sequence);
IO = IO{:}';
% Get analysis windows 
if(ischar(CONTROL(sf).stim_time))
   pspt = CONTROL(sf).stim_time;
   p = find( pspt~= ' ');
   pspt = pspt(p);
   stim_time=seq_parse(pspt);; % get the stimulus array - usually 2 times....
   stim_time = stim_time{:};
else
   stim_time = CONTROL(sf).stim_time;
end;

% get the psp definition array. For field potentials, we have the following
% sets of times: (layer 2 DCN pf FP)
% It is assumed that the same response is repeated with similar latencies for each stimulus...

if(ischar(CONTROL(sf).psp_time))
   pspt = CONTROL(sf).psp_time;
   p = find(pspt~= ' ');
   pspt = pspt(p);
   psp_time=seq_parse(pspt); % get psp definition array : for 
   psp_time = psp_time{:};
else
   psp_time = CONTROL(sf).psp_time;
end;

% Get sign information
if(ischar(CONTROL(sf).sign_list))
   pspt = CONTROL(sf).sign_list;
   p = find( pspt~= ' ');
   pspt = pspt(p);
   sign_list = seq_parse(pspt);; % get the stimulus array - usually 2 times....
   sign_list = sign_list{:};
else
   sign_list = CONTROL(sf).sign_list;
end;


ts1=number_arg(stim_time(1,1)); % get the baseline duration (most of time before stimulus)
ts2=number_arg(stim_time(1,2)); % get the baseline duration (most of time before stimulus)

tpb=number_arg(psp_time(1:2:end));
tpe=number_arg(psp_time(2:2:end));

% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
olddr = DFILE.rate(1);
k = find(DFILE.rate < 50);
DFILE.rate(k) = 100;
RATES = (DFILE.rate * DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec
RATES=ones(length(RATES), 1)*0.2;

tfac = DFILE.rate(1)/olddr;
time = time*tfac;

% calculate measurement times
fp_base = floor((ts1./RATES - 0.5)*0.75);
for i = 1:length(RATES)
   fpb(i,:) = floor(tpb./RATES(i));
   fpe(i,:) = floor(tpe./RATES(i));
end;
if(length(stim_time) > 1)
   dt = ts2-ts1;
   tpb2 = tpb + dt;
   tpe2 = tpe+dt;
   for i = 1:length(RATES)
      fpb2(i,:) = floor(tpb2./RATES(i));
      fpe2(i,:) = floor(tpe2./RATES(i));
   end;
end;

% measure holding current and "rmp"
for i = 1:records
   hold_cur(i) = mean(CURRENT(i,1:fp_base(i))); % save in array for later usage
   rmp_volt(i) = mean(VOLTAGE(i,1:fp_base(i)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(rmp_volt);

[vsmo, time] = FP_artsupp(VOLTAGE, DFILE, sf);

% but first smooth the voltage out a bit
for i = 1:records
   fsamp = 1000/RATES(i); % get sampling frequency
   fco = 1000;		% cutoff frequency in Hz
   wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
   if(wco < 1) % if wco is > 1 then this is not a filter!
      [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
      vsmo(i,:) = DigitalFilt(b, a, vsmo(i,:)')'; % filter all the traces...
   else
      vsmo(i,:) = vsmo(i,:);
   end
end

stdn = zeros(records, 1);
for i = 1:records
   base = mean(vsmo(i,1:fp_base(i)));
   stdn(i) = std(vsmo(i,[1:fp_base(i)]));
   vsmo(i,:) = vsmo(i,:) - base;
end;
stdn = sqrt(sum(stdn.^2)/records);
QueMessage('FP_TC - Measuring FP peak voltages and P3 slope for stim 1');

MAX_S = zeros(records, size(fpb,2));
MAX_S2 = MAX_S;
twinpulse = 0;
if(length(stim_time) > 1 & all(fpe2 <= size(vsmo,2)))
   twinpulse = 1;
end;

for i = 1:records
   for j = 1:size(fpb,2)
      if(sign_list(j) < 0)
         [npk, npi] = min(vsmo(i,fpb(i,j):fpe(i,j)));
         if(j == 1)
            [ppk, ppi] = max(vsmo(i,(fpb(i,j)+npi-1):fpe(i,j))); % neg to next positive (N1)
         else
            [ppk, ppi] = max(vsmo(i,fpb(i,j):(fpb(i,j)+npi-1))); % first pos to main neg. (N2)
         end;
         MAX_S(i,j) = abs(npk-ppk);
      else
         MAX_S(i,j) = max(vsmo(i,fpb(i,j):fpe(i,j))); % simply positive
      end;
      if(twinpulse)
         if(sign_list(j) < 0)
            [npk, npi] = min(vsmo(i,fpb2(i,j):fpe2(i,j)));
            if(j == 1)
               [ppk, ppi] = max(vsmo(i,(fpb2(i,j)+npi-1):fpe2(i,j))); % neg to next positive (N1)
            else
               [ppk, ppi] = max(vsmo(i,fpb2(i,j):(fpb2(i,j)+npi-1))); % first pos to main neg. (N2)
            end;
            MAX_S2(i,j) = abs(npk-ppk);
         else
            MAX_S2(i,j) = max(vsmo(i,fpb2(i,j):fpe2(i,j)));
         end;
      end;
      
      
   end;
end;
%   MAX_P3S(i) = max(diff(vsmo(i,fp_p3b(i):(fp_p3b(i)+pk)'))) / RATES(i); % max slope on rising phase
%end;
for i = 1:records
   a = find(MAX_S(i,:) < stdn);
   allsignal{i} = a;
end;

maxiter = 500;
maxio = max(IO);
minio = min(IO);
for j = 1:size(fpb,2)
   mxs = max(MAX_S(:,j));
   sl = (mxs-min(MAX_S(:,j)))/maxio;
   fit_bv(j,:) = [0:maxio];
   [fp2, chisq2, niter2, fit_pki2, err2, dep2] = mrqfit('sum_boltzmann', ...
      [min(MAX_S(:,j)) mxs*0.5 sl maxio/5 mxs*0.5 sl maxio*0.5], IO, MAX_S(:,j), [], ...
      [0 1 1 1 1 1 1], [-1, -mxs*0.5, 0.0, -maxio*0.5, 0, 0, minio*2], ...
      [-1, mxs*2, 10, maxio*1.5, mxs*2, 10, maxio*1.5], maxiter);
   fit_ba_double(j,:) = sum_boltzmann(fit_bv(j,:), fp2);
   fit_ba_dbl(j,:) = fp2;
   if(fp2(2) < fp2(5))
      fpx = fp2(1:4);
   else
      fpx = fp2([1,5,6,7]);
   end;
   [fb, chisqb, niterb, fit_b, errb, depb] = mrqfit('sum_boltzmann', [min(MAX_S(:,j)) mxs sl*2 maxio], ...
      IO, MAX_S(:,j), ...
      [], [1 1 1 1], [-1, 0, 0.0, -maxio*0.5], [1, mxs*2, 10, maxio*1.5], maxiter);
   fit_ba(j,:) = sum_boltzmann(fit_bv(j,:), fb);
   f1{j} = fb;
   f1d{j} = fp2;
   
   
end;

if(twinpulse)
   for j = 1:size(fpb,2)
      mxs2 = max(MAX_S2(:,j));
      sl = (mxs2-min(MAX_S2(:,j)))/maxio;
      fit_bv2(j,:) = [0:maxio];
      
      [fp2_double, chisq2, niter2, fit_pki2, err2, dep2] = mrqfit('sum_boltzmann', ...
         fit_ba_dbl(j,:), IO, MAX_S2(:,j), [], ...
         [1 1 1 1 1 1 1], [-1, -mxs2*0.5, 0.0, -maxio*0.5, 0, 0, minio*2], ...
         [1, mxs2*2, 10, maxio, mxs2*2, 10, maxio*1.5], maxiter);
      f2d{j} = fp2_double;
      fit_ba_double2(j,:)=sum_boltzmann(fit_bv2(j,:), fp2_double);
      
      if(fp2_double(2) < fp2_double(5))
         fpx = fp2_double(1:4);
      else
         fpx = fp2_double([1,5,6,7]);
      end;
      [fb, chisqb, niterb, fit_b, errb, depb] = mrqfit('sum_boltzmann', [min(MAX_S2(:,j)) mxs2 sl*2 maxio], IO, MAX_S2(:,j), ...
         [], [1 1 1 1], [-1, 0, 0.0, -maxio*0.5], [1, mxs2*2, 10, maxio*1.5], maxiter);
      fit_ba2(j,:) = sum_boltzmann(fit_bv2(j,:), fb);
      f2{j} = fb;
   end;
end;

FP.S1=MAX_S;
FP.S2=MAX_S2;
FP.IO=IO;
FP.S1_fitx=fit_bv;
FP.S1_fity=fit_ba;
FP.S1_fity2=fit_ba_double;
FP.S2_fitx=fit_bv2;
FP.S2_fity=fit_ba2;
FP.S2_fity2=fit_ba_double2;
FP.S1_par=f1;
FP.S1_par2=f1d;
FP.S2_par=f2;
FP.S2_par2=f2d;
FP.stdn = stdn;


CONTROL(sf).FP = FP; % copy the structure over

QueMessage('FP_IO - analysis complete');

na = size(MAX_S,2);
for i = 1:na
   IOX(:,i) = IO;
end;

% plot if figure is set
h = findobj('Tag', 'FP_IO'); % check for pre-existing window
if(isempty(h)) % if none, make one
   h = figure('Tag', 'FP_IO',...
      'Name', 'FP IO Analysis', ...
      'Units', 'normalized', ...
      'MenuBar', 'none', ...
      'WindowButtonMotionFcn', 'datac(''mouse_motion'', gcbf);', ...
      'WindowButtonDownFcn', 'datac(''mouse_down'', gcbf);', ...
      'WindowButtonUpFcn', 'datac(''mouse_up'', gcbf);', ...
      'NumberTitle', 'off');
end
figure(h); % otherwise, select it
clf; % always clear the window...
%Command Menu
uimenu('Label', 'Close &Window', 'Position', 1, 'callback', 'close(findobj(''Tag'', ''FP_IO''));' );
uimenu('Label', '&Print', 'Callback', 'print;');
uimenu('Label', 'Print&Setup', 'Callback', 'printdlg;');
fsize = 7;
msize = 3;
msg_pos = [0.85 0.90 0.15 0.07];
tmax = max(time);    

lab = {'N1', 'N2', 'P3', 'P4'};
npm = size(MAX_S, 2); % number of measurments
spc = 0.03;
wid = 1/npm - spc;
for i = 1:npm
   subplot('Position',[spc+(i-1)*(wid+spc),0.550,wid,0.30]);
   hold on;
   P = plot(IOX, MAX_S(:,i), 's', 'Markersize', msize);
   plot(fit_bv, fit_ba(i,:));
   plot(fit_bv, fit_ba_double(i,:), 'color', 'cyan');
   if(twinpulse)
      plot(IOX, MAX_S2(:,i), '+', 'Markersize', msize);
      plot(fit_bv, fit_ba2(i,:));
      plot(fit_bv, fit_ba_double2(i,:), 'color', 'magenta');
   end;
   
   set(gca, 'FontSize', fsize);
   ylabel('FP Amplitude (mV)');
   xlabel('Stim (uA)');
   grid;
   legend(P, char(lab(i)));
   u=get(gca, 'YLim');
   ux = get(gca, 'XLim');
   datac_setcrosshair(gca, 'FP_info', 'uA', 'mV', msg_pos);
   
end;

PP = zeros(records, size(MAX_S, 2));
for i = 1:records
   for j = 1:size(MAX_S, 2)
      if(MAX_S(i,j) < stdn)
         PP(i,j) = NaN;
      else
         PP(i,j) = MAX_S2(i,j)/MAX_S(i,j);
      end;
   end;
end;
subplot('Position', [0.55, 0.05, 0.43, 0.40]);
plot(IOX, PP, 'o', 'Markersize', msize);
hold on;
%plot(fit_bv', fit_ba2'./fit_ba');
plot([0:maxio], fit_ba_double2'./fit_ba_double');
u=get(gca, 'YLim');
set(gca, 'YLim', [0.5,u(2)]);
grid;

datac_setcrosshair(gca, 'FP_info', 'uA', 'Fac', msg_pos);
%	set(gca, 'YLim', [-1 1]);
%plot Max Slope voltage 
subplot('Position',[0.07,0.050,0.43,0.40]);
%   plot(IO(1:records),MAX_P3S(1:records), '-r^', 'Markersize', msize, 'Markerfacecolor', 'r') % data in red
plot(time', vsmo');
set(gca, 'FontSize', fsize);
v=get(gca, 'YLim');
vmax=max(abs(v));
%  set(gca, 'Ylim', [-vmax vmax]);
grid;
ylabel('mV)');
datac_setcrosshair(gca, 'FP_info', 'ms', 'mV', msg_pos);


% plot textual information abstracted from analysis...
subplot('Position',[0.07,0.85,0.73,0.13])
r_uncomp = CONTROL(sf).access;
axis([0,1,0,1])
axis('off')
text(0,0.75,sprintf('%-12s R[%d:%d] %-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8);
text(0,0.5,sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 7);
text(0,0.25,sprintf('Noise: %8.3f', stdn));

orient landscape

% control printing and closing of window for automatic runs
% f = 1 creates plot and leaves it up
% f = 2 creates plot but closes it when done
% f = 3 creates plot and prints it and then closes it
% control printing and closing of window for automatic runs
% plot_flag = 0 creates plot and leaves it up
% plot_flag = 1 creates plot and prints to the printer
% plot_flag = 2 creates plot and prints it and then closes it
% checking the pflag box in datac main window is plot_flag = 1;
if (plot_flag > 0)
   h = findobj('Tag', 'FP_IO');
   hg = findobj('Tag', 'FP_info');
   set(hg, 'Visible', 'off');
   print(h);
   set(hg, 'Visible', 'on');
end
if plot_flag == 2  
   close
end

